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SECTION  1 


SUMMARY 

Unlike  a  conventional  radar  signal,  an  ultra-wideband  (UWB)  radar  signal  contains  a 
wide  spectrum  of  frequency  components  that  are  subject  to  dispersive  effects.  The  radiated  and 
received  signal  shapes  will  be  altered  not  only  by  frequency-dependent  losses  in  the  antenna 
conductors,  but  also  by  wavelength-dependent  reflections  from  geometry  transitions  in  the 
antenna.  These  antenna-dependent  effects  must  be  minimized  and  characterized  before  a 
meaningful  target  analysis  can  be  performed.  The  antenna  effects  can  in  principle  be  treated 
computationally  using  time-domain  integral  equations,  which  offer  several  advantages  over 
alternative  approaches  using  either  the  frequency  domain  or  finite  differences. 

The  solution  of  dispersive  time-domain  integral  equations  for  practical  three-dimensional 
problems  has  been  severely  limited  by  the  additional  computational  complexity,  which  typically 
increases  by  a  factor  of  Nt  (number  of  time  steps)  over  the  corresponding  ideal  problem  with 
perfect  electrical  conductors.  The  primary  result  of  our  work  is  the  development  of  an  efficient 
solution  whose  computational  complexity  increases  by  a  factor  of  only  log(Nt)  over  the 
corresponding  ideal  problem.  This  efficient  solution  will  facilitate  the  use  computational 
techniques  for  the  design  of  UWB  antennas  and  for  the  analysis  of  transmitted  and  received 
signals. 


Section  2  summarizes  dispersive  effects  in  metal  conductors.  We  conclude  that  the 
conventional  low-frequency  model  for  electrical  conductivity  is  an  excellent  approximation  for 
frequencies  of  interest  for  UWB  radar,  and  that  there  is  no  experimental  or  theoretical  evidence 
of  other  effects  that  need  to  be  included  until  the  frequency  exceeds  100  GHz. 

Section  3  summarizes  the  existing  treatment  of  time-domain  integral  equations  for  three- 
dimensional  electromagnetics.  Most  of  this  work  treats  non-dispersive  problems,  i.e.,  perfect 
electrical  conductors  or  perfect  dielectrics.  Although  time-domain  integral  equations  have  been 
formulated  for  dispersive  media,  i.e.,  lossy  conductors  or  dielectrics,  only  limited  solutions  for 
simple  geometry  have  been  achieved  because  of  the  computational  challenges. 

In  Section  4  we  present  an  efficient  method  for  the  solution  of  dispersive  time-domain 
integral  equations  in  three  dimensions.  We  use  surface  impedance  boundary  conditions  to  avoid 
any  increase  in  the  number  of  scalar  integral  equations  that  need  to  be  solved.  The  resulting 
equations  still  contain  additional  time-convolution  integrals  that  increase  computational  time 
requirements  by  a  factor  of  Nt.  However,  we  have  developed  a  simple  recursion  relation  for  the 
convolution  so  that  the  computational  complexity  only  increases  by  a  factor  of  order  log(Nt)  over 
the  corresponding  ideal  problem. 

In  Section  5  we  describe  in  detail  the  numerical  implementation  of  the  efficient  solution 
for  a  two-dimensional  problem  involving  a  lossy  planar  transmission  line.  We  determine  the 
incident  fields  for  a  wave  that  is  fed  at  the  center  of  the  line,  and  then  reduce  the  time  and  space 


1 


2119  JR 


integrals  to  summations  that  can  be  performed  numerically.  The  resulting  equations  are  then 
solved  numerically  by  time-stepping  as  for  the  ideal  time-domain  integral  equation.. 

In  Section  6  we  validate  our  approach  by  comparing  the  computational  and  analytical 
solutions  for  the  two-dimensional  planar  transmission  line,  both  lossy  and  ideal.  We  find  good 
agreement  between  the  computational  and  analytic  results  in  both  the  time  and  frequency 
domains. 

Sections  7  and  8  contain  code  listings  and  references,  respectively. 
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SECTION  2 


ELECTROMAGNETIC  DISPERSION  IN  METALS 

Soon  after  the  discovery  of  the  electron,  Drude  [1]  proposed  a  theory  of  electrical 
conductivity  based  upon  the  assumption  that  there  is  a  gas  of  free  electrons  in  a  metal.  Electrical 
resistance  arises  from  the  collisions  of  accelerated  electrons  with  the  background  atoms.  Drude  s 
phenomenological  theory  gives  an  expression  for  the  low-frequency  conductivity 

ne2  e0co2p 

cf0  —  ’ 

my  y 

where  e  and  m  are  the  charge  and  mass  of  the  electron,  n  is  the  electron  density,  y  is  the  collision 
rate,  So  is  the  vacuum  permittivity,  and  cop  is  the  electron  plasma  frequency  given  by 


Modern  theories  of  electrical  conductivity  (Bardeen  [2]  and  Kittel  [3]  )  show  that  that  the 
electron  collisions  are  dominated  by  scattering  from  thermal  vibrations  of  the  metal  atom  lattice 
rather  than  from  individual  metal  atoms,  but  the  above  phenomenological  expression  remains 
valid. 


At  higher  frequency  one  must  include  electron  inertia  as  well  as  collisions.  The 
phenomenological  equation  for  the  electron  velocity  then  becomes 

mv  +  mvy  =  e  E  exp(-icot) . 


Then  one  obtains  a  complex  conductivity 

ne2 

o  = - = - — , 

m(y  -  ito)  (y  -  ico) 


and  a  corresponding  dielectric  function 


v  co(co  +  iy)J 

This  dielectric  is  a  function  of  frequency  and  thus  causes  dispersion  in  wave  propagation.  We 
have  searched  the  literature  extensively,  but  have  not  found  any  theoretical  or  experimental  basis 
for  other  effects  that  should  be  included  in  a  high-frequency  conductivity  model  for  metals. 


e(co)  =  s0  =  s0 
1(0 


A  comparison  of  key  conductivity  parameters  for  three  common  metals,  copper  (commercial 
annealed),  aluminum  (type  6061-T6),  and  stainless  steel  (type  304)  is  presented  in  Table  2.1. 

One  would  expect  high-frequency  conductivity  effects  to  become  significant  at  frequencies 
comparable  to  the  collision  frequency,  which  is  about  1013  Hz  (far  infra-red).  Indeed,  Stratton  [4] 
quotes  experimental  data  showing  that  the  low-frequency  conductivity  is  valid  up  to  about 
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1013  Hz.  This  frequency  is  a  factor  of  more  than  100  higher  than  the  tens  of  GHz  that  is 
envisioned  for  UWB  radar.  For  frequencies  above  the  plasma  frequency  at  about  3  x  1015  Hz 
(ultraviolet),  the  dielectric  has  a  positive  real  component  and  waves  will  propagate  rather  than 
attenuate.  Kittel  [3]  presents  experimental  data  that  confirm  the  ultraviolet  transmission  window 
for  thin  metal  films. 


Table  2-1 .  High-frequency  conductivity  parameters  for  common  metals. 


Metal 

Low-frequency 
conductivity  Go 
(mho/m) 

Conduction 
electron  density 
(m-3) 

Collision 
frequency  y/27t 
(Hz) 

Plasma 

frequency  cop/27i 
(Hz) 

Copper  (comm,  anld.) 

5.80  x  107 

8.45  x  1028 

6.53  x  1012 

2.61  x  1015 

Aluminum  (6061-T6) 

2.50  x  107 

1.81  x  1029 

3.24  x  10'3 

3.82  x  1015 

Stainless  steel  (304) 

1.40  x  106 

1.70  x  1029 

5.44  x  10'4 

3.70  x  1015 

Thus,  the  simpler  low-frequency  conductivity  model  is  valid  for  metals  and  frequencies 
of  interest  for  UWB  radar.  In  the  remainder  of  this  report  we  will  assume  that  the  metals  are 
good  but  nevertheless  lossy  conductors.  Physically,  this  means  ®  «  y  «  cop,  so  that  the 
complex  dielectric  function  is  well-approximated  by 

s(co)  =  s0(l -a/i(080),  ct/cos0»1. 

We  will  also  suppress  the  zero  subscript,  so  that  a  will  always  represent  the  low-frequency  (or 
“dc”)  conductivity. 
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SECTION  3 


TIME-DOMAIN  INTEGRAL  EQUATIONS 


In  1968  Bennet  and  Weeks  [5,6]  introduced  a  novel  approach  for  solving  the  time-domain 
integral  equation  for  the  scattering  of  electromagnetic  waves  from  perfectly  conducting,  two-  and 
three-dimensional  bodies.  Instead  of  resorting  to  the  matrix  solution  typical  of  frequency-domain 
solutions,  they  used  a  time-stepping  approach  similar  to  that  used  for  solving  initial  value 
problems.  Early  reviews  by  Poggio  and  Miller  [7]  and  Mittra  [8]  provide  detailed  derivations  of 
the  basic  equations  and  extensive  details  of  the  computational  implementations.  A  more  recent 
review  by  Miller  [9]  summarizes  applications  of  the  time-domain  approach,  all  of  which  are 
limited  to  perfect  conductors.  Although  Tesche  [10]  investigated  the  time-domain  integral 
equation  for  lossy  conductors  in  1991,  the  first  numerically  rigorous  solution  for  a  three- 
dimensional  problem  was  reported  ten  years  later  by  Bluck  et  al  [11]. 


We  first  review  the  three-dimensional  time-domain  integral  equations  for  propagation  in 
a  linear,  isotropic,  homogeneous  medium.  Poggio  and  Miller  [7,  equ.  4.35  &  4.37]  give  the 
electric  and  magnetic  field  integral  equations  (EFIE  and  MFIE  respectively)  for  a  source-free 
volume  bounded  by  a  closed  surface  S  (part  of  which  may  be  at  infinity) 


_n 

471 


E  =  E" 


■  —  cfds'  fdt' 
4ttT  J 


R 


—  H  =  Hinc  + 
471 


—  cfds'  fdt' 
4ti  T  J 


R 


[4ln'*H'1- 

[I+IAlfc 

|_R  c  at'JL 

1  1  d\ 
R  +  c  dt'J 

Here  bold  type  denotes  vector  quantities,  and 

E  =  E(r,t)  is  the  electric  field  at  the  observation  point; 

E’  =  E(r’,t')  is  the  electric  field  at  the  source  point  on  the  surface  S; 

Einc  =  Einc(r,t)  is  the  incident  electric  field  generated  by  sources  outside  of  S,  typically  a 
plane  wave  from  infinity  or  a  current  at  an  antenna  feed; 

Similar  conventions  hold  for  H,  H',  and  H,nc; 

Q  =  Q(r)  is  the  solid  angle  subtended  by  the  medium  at  the  observation  point,  e.g.  271  on 
a  smooth  surface  S  or  4tc  in  the  enclosed  volume; 

8  is  the  Dirac  delta  function; 

8  and  p  are  the  electric  permittivity  and  magnetic  permeability,  respectively; 
c  =  (ep)’05  is  the  speed  of  light  in  the  medium; 

R  =  r  -  r1,  R  =  |R|,  and  R  is  the  unit  vector; 
n1  is  a  unit  normal  to  S  pointing  into  the  enclosed  volume; 
is  a  principal  value  integral  on  the  surface  S. 

Note  that  the  time  integral  is  trivial  for  non-dispersive  media,  with  the  source  quantities  simply 
evaluated  at  the  retarded  time  f  =  t  -  R/c. 
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When  the  surface  S  coincides  with  a  perfect  electrical  conductor  (PEC),  the  electric  and 
magnetic  fields  must  satisfy  boundary  conditions  given  by 

nxE(r,t)  =  0  and  n-H(r,t)  =  0. 

For  perfect  electrical  conductors  it  is  also  useful  to  define  the  surface  current  and  charge  density 

Js(r,t)  =  nxH(r,t) 

ps(r,t)  =  sn  -E(r,t) , 


which  satisfy  a  conservation  relation 

4-ps(r,t)  +  VJs(r,t)  =  0. 

at 


Taking  the  cross-product  with  the  normal  n  to  the  surface  S,  and  using  Q  —  2n  for  r  on  S?  the 
EFIE  and  MFIE  can  then  be  written  as 


nxE-^D^^ds'j^Js^t')- 


j_  _LA 

R2  +  Rc  8t' 


PsC1*'^') 


I  t'=t— R/c 


Js(r,t)  =  2nxHinc(r,t)  +  ^-x<j[ds' 


ffi  |  i  a 

{|_R2  +Rc  dt' 


Js 


t'=t-R/c 


Either  one  of  these  equations  suffices  to  determine  the  two  scalar  components  of  Js,  which  in 
turn  determines  the  electric  and  magnetic  fields  in  the  volume  enclosed  by  S.  Miller  notes  that 
the  MFIE  is  best  suited  for  smooth  closed  conductors,  while  the  EFIE  is  applicable  to  more 
general  bodies  including  thin  wires  and  shells. 


When  the  surface  S  coincides  with  a  lossy  conductor,  then  one  must  solve  a  second  pair 
of  integral  equations  for  the  lossy  medium,  and  then  match  boundary  conditions  across  S.  Bluck 
et  al  [11]  write  the  EFIE  and  MFIE  for  a  lossy  conductor  (using  the  same  convention  for 
observation  and  source  points)  as 


— E  =  — cf  ds'  fdt'J 
471  471  *  J 


pn  x 


5HF 

at' 


iG  -  (n'  x  E')  x  V'G  -  (n'  •  E" ) x  V'G 


_Q 

471 


H=  —  cf  ds'  fdt'J — 
471*  J 


,  8E! 

on xE  +snx — - 

at' 


G-(n'xH')xV'G-(n'-H')xV'G  , 


G(r,t;r',t')  =  exp(-px) 


S(t-R/c)  ]  P/c  t  0(x -R/c) 


R 


■R  /c 


Here  G(r,t;r',t')  is  the  Green’s  function  given  by  Morse  and  Feshbach  [12],  the  material  properties 
now  refer  to  the  lossy  medium,  and 

p  =  g/2s  is  the  characteristic  decay  rate; 

n1  is  a  unit  surface  normal  directed  away  from  the  lossy  body; 

t  =  t  —  t'; 
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Ij  is  the  first-order  modified  Bessel  function  of  the  first  kind; 

0  is  the  Heaviside  unit  step  function. 

The  boundary  conditions  at  the  interface  S  between  the  ideal  dielectric  (subscript  1)  and  the  lossy 
conductor  (subscript  2)  are  given  by 

s,  n-E,(t)  =  e2n-E2(t)  +  a2  J^du  n-E2(u) 

n  (B,  -B2)  =  0 


nx  (E,  -E2)  =  0 


nx(H,  -H2)  =  0. 

Bluck  et  al  [11]  have  computed  the  scattering  of  a  gaussian  plane  wave  by  a  lossy  dielectric 
sphere,  and  have  obtained  good  agreement  with  the  analytical  solution 

Nevertheless,  the  computational  complexity  of  the  lossy  time-domain  integral  equation 
has  prevented  its  application  to  practical  problems.  This  complexity  arises  primarily  for  three 
reasons.  First  of  all,  the  lossy  problem  requires  a  true  time  integration  (rather  than  a  delta- 
function)  at  each  space  point  for  each  time  step,  which  increases  the  computational  complexity 
by  a  factor  of  order  Nt,  the  number  of  time  steps.  Second,  the  lossy  problem  requires  the 
computation  of  12  field  components  (3  electric  and  3  magnetic  for  both  media)  at  each  point  in 
space  and  time,  compared  to  only  two  components  of  Js  for  the  PEC  problem.  Finally,  the  lossy 
Green’s  function  requires  the  evaluation  of  exponential  and  Bessel  functions,  rather  than  the 
simpler  1/R  of  the  PEC  problem.  In  the  following  section  we  propose  an  alternative  approach 
that  greatly  mitigates  the  complexity  introduced  for  the  first  two  reasons. 
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SECTION  4 


EFFICIENT  SOLUTION  FOR  LOSSY  CONDUCTORS 


In  this  section  we  propose  an  efficient  solution  for  the  dispersive  time-domain  integral 
equation  that  applies  when  the  conductivity  is  high  enough  that  the  electromagnetic  skin  depth  is 
much  smaller  than  any  geometric  scale  length.  For  small  skin  depth  the  problem  can  be  greatly 
simplified  by  the  use  of  an  impedance  boundary  condition  at  the  surface  of  the  material.  We  then 
use  an  improved  exponential  approximation  for  the  surface  impedance  and  develop  a  recursion 
relation  that  greatly  reduces  the  computation  involved  in  performing  the  convolution  integral. 

The  surface  impedance  boundary  condition  (SIBC)  is  attributed  to  Leontovich  [13],  but 
Senior  [  1 4]  provides  a  more  accessible  exposition.  Mitzner  [  1 5]  first  used  the  SIBC  in  a 
frequency-domain  integral  equation  approach  to  scattering  from  a  body  of  finite  conductivity. 
Tesche  [10]  first  applied  the  SIBC  to  a  time-domain  integral  equation,  but  his  non-rigorous 
approach  resulted  in  errors  that  we  have  corrected  below. 


We  assume  the  conductor  is  a  good  conductor,  i.e.,  c»coe  for  all  frequencies  of  interest, 
and  that  the  skin  depth 

\  jxcoa 


is  negligible  compared  to  the  wavelengths  of  interest,  to  the  radii  of  curvature,  and  to  any  other 
geometric  scale  lengths  in  the  problem.  We  have  rewritten  Senior’s  [14]  frequency-domain 
SIBC  in  terms  of  a  normalized  impedance  Z  that  is  a  function  of  the  Laplace  variable  s, 


E  -  (n  •  E)n 


ZnxH 


Here  n  is  the  unit  surface  normal  pointing  away  from  the  conductor,  and  the  subscript  c  refers  to 
conductor  properties.  Taking  the  cross-product  with  n,  the  SIBC  can  be  rewritten  as  a  surface 
Ohm’s  law 


To  obtain  a  boundary  condition  for  the  normal  component  of  the  magnetic  field  we 
rewrite  Senior’s  result  using  the  normal  derivative  of  H  as 


Y  =  - ,  P  =  — , 

y  s(s  +  2p)  H  2s  c 

where  Y  is  a  normalized  pseudo-impedance.  We  then  use  V  •  H  =  0  and  V  •  n  =  Vn  =  0  for 
Cartesian  coordinates  to  express  the  normal  derivative  in  terms  of  the  surface  current 

V  x  (n  x  H)  =  H  Vn  -  n  VH  +  nV  H  -  HV  n  =  -  n  VH  . 

Thus,  we  obtain  desired  boundary  condition  on  the  normal  component  of  H, 


n  H  =  -  c  J— —  Y  n  •  V  x  Js . 

rJ 


Since  the  above  SIBCs  on  n  x  E  and  n  •  H  are  in  the  frequency  domain,  we  must  perform 
the  inverse  Laplace  transforms  to  obtain  the  corresponding  boundary  conditions  in  the  time 
domain.  Using  standard  Laplace  transform  tables  we  obtain 

Z(t)  -  5(t)  -  p  e-pi  [I0  (Pt)  -  I0(pt)]  ^  6(1)  -  PZ,  (t) 


Y(t)  =e“pt  I0(Pt), 


where  Io  and  L  are  modified  Bessel  functions  of  the  first  kind.  The  time-domain  boundary 
conditions  then  involve  convolution  integrals 


nxE 


s 


n  •  H  =  -  c  j— —  Y  *  n  •  V  x  Js 

V  8c  P 

where  f  *  g  =  dt  f  (u)  g(t  -  u) . 

In  the  limit  of  perfect  conductivity  ( p  — »  oo )  these  boundary  conditions  reduce  to  those  for  a 
PEC,  i.e.,  nxE  =  n-  H  =  0. 


To  obtain  the  time-domain  integral  equations  for  a  region  bounded  by  a  lossy  conductor, 
we  first  take  the  cross-product  with  n  of  the  first  EFIE  and  MFIE  of  the  preceding  section,  giving 


JQ 

471 

Q 

— r 
4tt 


fds'-j 

'i  i  a" 

[[n'  x  E'l  x  R  +  [n'  •  E']r]1 

Js  r  j 

at' 

_r  c  at'. 

LL  J  J 

cfds'  — < 

U[n'xE']  + 

'i  i  a' 

|n'  x  H']  x  R  +  [n'  •  H']r]1 

i  r 

[  at' 

_r  c  at'. 

J 

t'=t-R/c 


t'=t-R/c 


Next  we  insert  the  above  time-domain  boundary  conditions  and  use  the  charge  conservation 
relation  between  Js  and  ps  to  obtain 
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Since  either  equation  suffices  to  determine  the  surface  current  density  Js,  we  have  reduced  the 
number  of  scalar  field  quantities  from  twelve  to  two.  Thus,  we  have  eliminated  the  second 
complexity  issue  raised  at  the  end  of  the  preceding  section.  Furthermore,  we  have  corrected  the 
error  in  Tesche’s  result,  which  omitted  the  convolution  Z'  *  Jg  term  for  source  points  on  S. 

These  integral  equations  still  involve  a  true  time  integration  in  the  convolution  term.  The 
computational  complexity  can  be  greatly  reduced  if  the  convolution  integral  can  be  recursively 
updated.  In  1990  Luebbers  et  al  [16]  demonstrated  that  a  recursive  convolution  was  possible  for 
a  time-domain  dielectric  susceptibility  given  by  a  simple  exponential  function.  The  recursive 
convolution  integral  was  soon  applied  to  a  finite-difference  time-domain  SIBC  by  Maloney  and 
Smith  [17],  who  used  Prony’s  method  to  approximate  the  surface  impedance  by  a  sum  of 
exponentials.  Oh  and  Schutt-Aine  [18]  fitted  a  minimax  rational  approximation  to  the 
frequency-domain  impedance,  then  decomposed  the  fit  into  a  sum  of  simple  poles  and  obtained 
an  exponential  sum  representation  for  the  surface  impedance  in  the  time  domain.  Unfortunately, 
Prony’s  method  requires  a  large  number  of  exponential  terms  for  a  good  approximation  over  a 
large  number  of  time  steps,  while  Oh’s  approach  gives  a  poor  approximation  in  the  time  domain. 
These  problems  arise  primarily  from  the  huge  variation  in  the  time  domain  impedance  Zi(u), 
which  asymptotically  approaches  u'3/2  for  large  arguments  u. 

To  obtain  an  improved  exponential-sum  approximation  using  a  limited  number  of  terms, 
we  have  used  a  different  approach  to  determine  the  exponents  and  coefficients  in  the  fit 

M 

Zl(Pt)  =  SC">eXP(-®mPt)' 

m=] 

Specifically,  we  have  chosen  the  exponents  (Oj  to  span  the  range  [(3N,PA)‘',  3  PA]  in  a 
logarithmically  uniform  manner,  where  A  is  the  time  step  and  Nt  is  the  number  of  steps  in  the 
numerical  implementation.  Then  the  coefficients  Cj  are  determined  by  a  least  squares  fit  to 
Zj(nPA)  over  the  range  1  <  n  <  N(,  again  using  logarithmically  uniform  spacing  for  n.  As  shown 
in  Figure  4-1,  this  approach  provides  an  excellent  minimax  approximation  to  Z\  over  the  entire 
time  range  for  the  computation.  The  number  of  terms  required  depends  upon  the  desired  relative 
error  (8Z/Z)  as  well  as  the  number  of  time  steps,  and  is  approximately  given  by 

M»21ogjN,  A); 

M  =  1 0  terms  give  a  relative  error  of  <0.2%  over  1 03  time  steps. 
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Figure  4-1.  Relative  error  for  a  10-term  exponential  fit  to  the  normalized  impedance  Zi(n(3A). 

Before  using  the  exponential  fit  in  the  convolution  integral,  we  have  found  it  helpful  to 
separate  out  the  time  interval  [0,A]  and  to  rewrite  the  scalar  convolution  integral  as  the  sum 
(suppressing  the  spatial  argument  for  simplicity) 

E(t)  s  Z  *  J  =  du  Z(t  -  u)  J(u)  +  j*_A  du  Z(t  -  u)  J(u) 

=  <f_A  du  Z(t  -  u)  J(u)  -  f A  du  p  Z,  (P(t  -  u))  J(u) 

=  £dv[5(v)-pZ,(pt)]  J(t-v)-]Tcm  £  du  p  exp(-comP(t  -  u))  J(u) 

m=l 

M 

-E.M-JXO)  . 

m-1 

Assuming  that  the  time  step  A  is  chosen  small  enough  that  J(t)  is  nearly  constant,  we  can  pull  J(t) 
outside  the  Eo(t)  integral  and  use  [19,  equ.  6.61 1.4] 

jJ°duZ,(u)  =  l 

to  obtain  an  accurate  small  difference  between  two  large  terms 

E0  (t)  =  J(t)  -  J(t)  J^A  dv  Z,  (v)  =J(t)|°AdvZ,(v)  . 
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We  then  apply  the  good  conductor  approximation  pA  »  1  use  asymptotic  approximations  for  the 
modified  Bessel  functions 


I0(x) : 

I,(x) 


e* 


,19  75 

1  H - 1 - —  + 


V2 


7TX 


8x  128x2  1024x3 

15  105 


+  .... 


>-A_ 


8x  128x  1024x- 


for  x»l 


for  x»l 


to  obtain  a  simplified  result  for  the  initial  tern  in  the  convolution 


E0(, )*!(<)  £dvZl(v)*JJ=[l  +  A 


8pA  128(PA) 


■  +  ... 


For  the  remaining  convolution  terms  Ei(t)  through  EM(t)  we  can  derive  a  recurrence 
relation.  Taking  J(t)  to  be  nearly  constant  within  each  time  step,  we  write  for  n  >  1 

Em(nA)  s  ')Adu  pcm  exp(-p(nA-u))  J(u) 

=  cm  ')Adu  p  exp(-comP(nA - u))  J(u) 

=  cmZ  JTnAdu  P  exp(-comP(nA-u))  J(kA) 

=  —  yexp(-©mp(n-k)A)[l-exp((-compA)]  J(kA) 

©m  tt 

=  exp(-compA)-^-y  exp(-comP(n  -k)A)  [1  -exp((-compA)]  J(kA) 

©m  tt 

+  -^-exp(-©„,pA)[l  -  exp((-(ompA)]  J((n  -1)A) 

©m 

=  exp(-co  PA)Em  ((n  -  1)A)  +  -^exp(-compA)[l  -  exp((-conipA)]  J((n  -  1)A) 

©m 


The  recurrence  relation  involves  only  the  values  of  the  convolution  term  Em  and  the  surface 
current  J  at  the  previous  time  step  (n-l)A.  A  similar  recurrence  relation  can  be  developed  for  the 
convolution  of  Y(t)  with  the  current.  Since  there  are  only  M  terms  that  need  to  be  updated  at 
each  space-time  point,  we  have  reduced  the  added  computational  complexity  of  the  convolution 
from  a  factor  of  order  Nt  to  a  factor  of  order  log(Nt)- 
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SECTION  5 


IMPLEMENTATION  FOR  A  LOSSY 
2-D  PLANAR  TRANSMISSION  LINE 

In  this  section  we  illustrate  the  efficient  solution  for  a  specific  problem  and  develop  the 
numerical  equations  in  detail.  We  have  chosen  a  two-dimensional  (2D)  rather  than  a  three- 
dimensional  problem,  not  only  because  of  our  limited  computational  resources,  but  also  because 
we  desire  a  problem  with  an  analytical  solution  that  can  validate  the  efficient  MFIE  solution. 

The  lossy  2D  planar  transmission  line  illustrated  in  Figure  5-1  has  a  well-established 
analytical  solution  [20].  The  line  consists  of  a  free-space  region  located  between  two  thick 
planes  of  a  lossy  conductor  with  dielectric  constant  £o  and  conductivity  a.  The  volumes  are 
separated  by  a  gap  a  in  the  x-direction,  the  y-coordinate  ignorable,  i.e.,  did y  =  0 ,  and 
propagation  is  in  the  z-direction.  All  fields  are  zero  before  t  =  0,  when  a  driving  current  J 0  (t)x 
is  applied  across  the  gap  at  the  location  z  =  0.  Waves  will  propagate  in  the  +  z  directions  with 
surface  current  ±  Js  (z,t)z  on  the  upper  surface  at  x  =  +a/2  and  +  Js(z,  t)z  on  the  lower  surface  at 

a  =  -a/2.  Because  of  this  symmetry  we  only  need  to  solve  Js(z,t)  for  observation  points  at  x  =  a/2 
and  z  >  0.  The  actual  field  components  are  Hy,  Ex,  and  Ez ,  but  these  are  secondary  variables  that 
are  derivable  from  the  solution  Js(z,t).  We  numerically  implement  only  the  MFIE,  which  is 
preferable  to  the  EFIE  for  this  type  of  problem.  Taking  ec  =  s  =  So,  pc  =  B  =  Ho  ,  =  2n,  and 

V  x  zJ(z,t)  =  0 ,  the  MFIE  for  the  2D  planar  problem  becomes 


z  =  0 


Figure  5-1 .  Geometry  for  2D  planar  transmission  line. 


The  incident  field  Hinc(z,t)  results  from  the  current  applied  at  the  plane  z  =  0  in  the 
absence  of  the  conductors  at  x  =  ±a/2.  Although  we  could  have  generated  the  incident  magnetic 
field  by  a  current  applied  on  a  conductor  closing  the  gap  at  z  =  0,  it  is  numerically  easier  to  apply 
twice  the  current  across  free  space,  which  gives  the  symmetric  waves  traveling  ±  z  directions. 
The  incident  magnetic  field  Hmc(r,t)  is  given  by  the  vector  potential 

Hinc (r,t)  =  Vx  Ainc(r,t) 

where  the  time-domain  vector  potential  for  a  sheet  current  2 Jo  at  z1  =  0  and  |xf|  <  a/2  is  given  by 


^  a/2  oo  oo 

Ainc  (r,t)  =  —  Jdx'Jdy'Jdt 

^-a/2  -oo  0 


t 


5(t'-t  +  R/c) 
R 


2J0(t)x  . 


The  incident  surface  current  on  the  surface  at  x  -  a/2  is  then  given  by 


a/2 


J inc (z, t)  =  n  x  H inc  (r , t)  =  — —  x  fdx'  fdy' fdt' V 

4n  J  J 


-a/ 2  -oo  0 


R 


a  a/2  oo  oo 

=-^x  K'M®' 


n  -a.  o 

a  /2  oo  oo 


8(t'-t  +  R/c)  t  5'(t'-t  +  R/c) 


=  |-x  Jdx'Jdy'Jdt' 

-a/2  -oo  0 
a  a  /  2  oo  oo 

KM*' 

Z,;t-a/2  -oo  0 

=  z  Jinc(z,t)  , 


R2 

'jpco ,  gj0(tf) 

R2  cR  at' 

Jo(t')  dJ0(t') 
R2  cR  at' 


cR 


x  x2J0(t') 
Rx  xJ0(t') 


(a/2-x')x  +  (y-y')y +  zz 
R 


—  5(t'-t  +  R/c) 
R 


xx5(t'-t  +  R/c) 


where  we  have  used  the  fact  that  that  the  (y-y1)  component  of  R  is  odd  and  vanishes  when  the  y 
integration  is  carried  out. 


Returning  to  the  vector  MFIE,  we  perform  the  vector  algebra  to  reduce  the  MFIE  to  a 
scalar  equation.  Since  the  unit  normal  n  points  into  the  gap,  the  first  vector  product  gives 

x'  =  +a/2 :  nx(J'xR)  =  (-x)x(zJ'x(y-y'  +  z-z'))  =  xxxJ'x(y-y')J'  =  0 
x'  =  -a/2 :  nx(J'xR)  =  (-x)x(-zJ')x(ax  +  y-y'  +  z-z'))  =  xxyaJ'  =  zaJ'  . 


The  second  vector  product  gives 

x'  =  +a  /  2 :  n  x  (n'  x  J')  =  (-x) x  (-x x  zJ')  =  (-x) x  (yj')  =  -zJ' 
x'  =  -a  /  2 :  n  x  (n'  x  J')  =  (-x)  x  (x  x  (-zJ'))  =  (-x)  x  (yJ')  =  -zJ' 

Thus  all  terms  have  only  z  components  and  we  obtain  a  scalar  integral  equation 

J(z,  t)  =  2 Jinc  (z,  t)  +  J ,  (z,  t)  -  J  2  (z,  t)  -  J3  (z,  t) 


where  the  individual  terms  are  given  by 
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<•  a/2  oo  oo 

Jinc(z,t)  =  —  fax'  fdy' fdt'-8(t'-t  +  R/c) 

271 -a/2  -t  0  R 


Jo(t') ,  gjp(t') 
R2  cR  dt' 


R2  =(|-x')2  +(y-y')2+z2 


_  i 


..00  00 

j  I  (Z,  t)  =  Jdz'  Jdy'jdt 


2tx 

if- 


—00  —00 
00  00 


R  +  cdt' 


7\2 


-J'(z',t')  ,R2  =  (-)2  +(y-y')2  +(z-z') 
R  2 


dt'^  -  -  JV~— ■{  Z'(t')  *  J'(z',(t')  L  R2  =  £)2  +  (y  -  y')2  +  (z -  z')2^ 


i  +  i  a  7V. 


J3(z,t)-  —  Jdz'  Jdy'jdt 


_f\  2 


—  00  —00 


R 


i^Z'(t')  *  J'(z',t')  ,  R2  =  (y  -  y')2  +  (z - z') 


[  c  dt 


Each  of  these  terms  contains  a  y’  integration  where  y'  appears  only  in  R.  Without  loss  of 
generality  we  take  y  =  0  and  change  the  integration  variable  from  y'  to  R 


,  8(t'-t  +  R/c) 
Rn 


dR 


8(t'-t  +  R/c) 


R"'1  ^R2  -(x-x')2  -  (z  -  z')2 


[c  1 1  - 1'  I]"'1  ^c2(t-t')2  -(x-  x')2  -  (z  -  z') 


=  KdR8(R-c(t-t')) 

/\  2  «b 


if  c(t - 1')  >  -x')2  +  (z- z')2 


rc|t-t'  |]—  ^  -  if  -  (x  -  o’  -  <z  -  z')J 

[0  otherwise  . 

This  result  is  consistent  with  the  2D  Green’s  function,  whose  time  dependence  has  an  inverse 
square  root  rather  than  a  delta  function  [12]. 

Thus  the  Jinc(z,t)  term  becomes,  after  performing  the  y1  integration, 


a/2  I 


Jinc  (r,  t)  =  —  Jdx'Jdt' 


Jo(t')  ,  1  ajp(t') 


2c 


x  >  xmin  =  maxi 


277xl  oJ  L  c2(t-t')2  c2(t-t')  a'  J  Vc2(t-t')2-(a/2-x')2-z2 

where  the  integral  limits  are  determined  by  the  requirement  that  c2(t-t')  2-(a/2-x‘)  -z  >0,  i.e., 

J  a/2-A/c2(t-t')2  -z2 
[  -  a/2 

This  result  can  be  further  simplified  by  changing  variables  to  u  =  c(t-t')  and  v  —  x  -a/2, 

z  r  i  Jn(t  —  u/c)  8J0(t  — u/c) 

J'nc(r,t)  =  —  fdu  oV  2  0  — 

no 

vmin 

F(u)  =  Jdx' 


Vu2  -  v2  -z2 


=  sin 


u  da 

J 

1  9 

(  a  1 

vmax  s  minj 

v 

V  max  / 

5 

}  -n/u2 
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Integrating  the  derivative  term  by  parts, 


Cl 

Jdn 


5J0(t-u/c) 


du 

=  -— J0(t-z/c)  +  fdu  j0(t-„/c)M-  fdu  J0(t-u/c)^^ 
2z  J  u  i  u  du 


F(u) 


u 


J0(t-u/c) 


F(u) 


net  ct 


U 


Jdu  J0(t-u/c) 


-Jz  z 


aF(u)  F(u) 

U  0U  U2 


For  u2  <  a2  +  z2  we  have  vmax  =  a  and  F  =  constant.  Thus  5F/du  is  non-vanishing  only  for 
u2  >  a2  +z2,  where 


5F(u)  d  . 
— —  = — sin 

au  au 


Vu2  -z2  ) 


-au 


au 


->/l -  a2  /(u2  -z2)  7(u2  -z2)3  (u2  -z2)Vu2  -z2 


Combining  the  above  results  we  obtain  the  incident  field 


1  Ct 
Jln‘(M)  =  -J„(<-z/c)-—  fdu 

)  TT  J 


Jo(t-u/c) 


*1  (u2-z2)Vu2-z2-a2 


a  =  Va2  +z2 . 


Thus  Jmc(z,t)  vanishes  for  z  >  ct,  because  Jo(t)  vanished  for  t  <  0.  Taking  Jo(t)  to  be  piecewise 
linear,  we  can  rewrite  the  incident  field  as  a  sum  of  separate  integrals  over  each  time  step 

Jinc(z,t)  =  ^-J0(t-z/c) — 2  [J0(t-u/c)-J0(t-u/c  +  A)][F,(u)-F,(u-cA)] 

2  71CA  U=a+CA 

- —  t  [(cA-u)J0(t-u/c)  +  u  J0(t-u/c  +  A)][F0(u)-F0(u-cA)] 

2ticA 


k  u=a+cA 


2zadu 


u 

F°(u)s  f  2  2  n — 2 — : 

(u2  -z2)vu2  —  z2  —  a'1 

U 

f  audu 

F,(u).  j 


■  =  sin 


-i 


^z2  -  a2  -2a2z2  /(u2  -z2)^ 
z2  +a2 


(u2  -z2)a/u2  -z2  -a2 


=  tan' 


u2-z2 


-1 


The  Ji(z,t)  term  becomes,  after  performing  the  y1  integration, 


zmax  t-a, 

J,(z,t)  =  -  fdz-  fdf 

TT  J 


zmin  0 


JCz'.Q 

f'A2 


1  5J(z',t') 


c 2 (t — t* ) 2  c2(t-f)  at’ 


V°2(t - 1')2  -a? 


The  finite  limits  on  the  z1  integration  arise  because  the  incident  current  Jinc(z,t)  and  consequently 
the  total  surface  current  J(z,t)  are  identically  zero  unless  t  >  z.  Physically,  the  sum  of  the 
propagation  time  of  the  incident  wave  from  the  origin  to  the  source  point  z1  plus  the  additional 
propagation  time  from  the  source  point  to  the  observation  point  z  must  be  less  than  t  in  order  for 
get  a  non-zero  contribution  to  J(z,t),  i.e., 
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It  is  then  straightforward  to  show  that  the  z1  integrand  vanishes  outside  the  interval  [zmjn,zmax], 
where 


7  = 

max 


^Tnin 


t2-z2-(x-x')2 
2(t  -  z) 

t2-z2-(x-x')2 
2(t  +  z) 


Next  we  change  the  time  variable  to  u  s  c(t-t')  and  integrate  by  parts  to  remove  the  (t-t1)' 
singularity, 


j,(Z,t)=-  rmaxdz'  r  a,dt'U 
*  to  Lc 

71  bmin  •ca,  L  u 


J(z',t') 


1  <3J(z',t') 


c2(t-t')2  c2(t-t')  at'  j  /c2(t_t,)2  _a2 


/c)  1  dJ(z',t-u/c)  1 


cu  at' 


u  -a, 


3  j^ma 

n  kmin 


I u2-ot2  fct  aj(z',t-u)  Vu2_a?  1 

2  J(z,t  — u)  +  I  du  +  - - - 

afu  *fca,  eft  ua,  u,/u2-aj 


dz’  r  du 

•tai  c^t  a2  Ju2  -af 


Next  we  take  J(z',t')  to  be  piecewise  linear  from  one  time  point  to  the  next  so  that  the  slope  is 
constant  within  each  time  interval  [t',t'+A],  Then  we  can  write  the  u-integral  as  a  sum  of 
integrals  with  aj  /  at  outside  of  each  time  integral, 


J,(z,t)  =  -  fmaxdz'  I 

71  JK  ■  - 


t-A  i/V 


J(z', t  —  u / c)  —  J(z', t  —  u / c  —  A)  fu  +  cA  u 


r  +  CA  u 

a,2-Ju2  -  a? 


_  3_  finax  t  A  J(zy,t-u/c)- J(z;?t-(u/c  + A))  Vu2 
7i  -  c  A  a? 


u  +  cA 


_  a  j^n 

7i 


dz'  I 


J(z',t-u/c)-J(z',t-(u/c  + A))  -y(u  +  cA)2  -a2  -^/u2  -a,2 


The  z'  integrand  is  a  complicated  function  of  z1  because  ai  depends  upon  z',  so  we  will  evaluate  it 
numerically  rather  than  analytically.  Using  the  trapezoidal  rule  with  a  step  size  cA  that  meets  the 
stability  criterion  8R  <  c8t ,  we  obtain 
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-max  t-A  ,  J(u  +  cA)2  -a?  -Ju2  -a 

J,(z,t)  =  -  I  <(»(z')  I  [j (z', t  -  u / c)  -  J (z', t  -  (u / c  +  A))]  *  - - 


z 

&  max 


nz'=z-  u  =  a, 

min  i 


a: 


where  we  have  omitted  the  0.5  factor  for  the  end  terms  because  the  integrand  actually  vanishes 
there.  We  have  also  included  the  factor 


<Kz')  = 


1  for  z'  >  0 
0  for  z'  =  0 
- 1  for  z'  <  0 


since  J(z’,tl)  changes  sign  at  z1  —  0,  which  also  results  in  a  null  contribution  at  z1  —  0. 

The  J2(z,t)  term  becomes,  after  performing  the  y1  integration, , 

i  zmax  t-al  * 

j2(z,t)=-  jdz'  J(z%t-> . 

"zrnin  0  C5t  VC2(t-t  )2  -a, 2 

where  a,  =J(x-  x')2  +  (z  -  z')2  . 

Then  writing  E(z',t’)  =  Z(t')*  J(z',t')  to  simplify  the  notation,  the  evaluation  of  J2(z,t)  proceeds  in 
a  manner  similar  to  Ji(z,t)  except  that  no  integration  by  parts  is  needed, 

,  (*“«,  <5E(z',t')  c 


,!(z’t)4r"“dz'i  a’dt  ca 

*^min 

=  1  fn,axdz'  f*  duaE(Z'’t_U/c)  1 
7i  .ba, 


-x/c2(t - 1')2  -a2 


•^/u2  -  a 


= i  fmax dz' c(tfA)  E(z,»t-u/e)-E(z,»t~u/c~A)  r +  du _ L= 

11  E.  Z„,cc,  cA  *■ 

■if 

71  K 


:  c(t-A)  E(z',t-u/c)-E(z',t-u/c-A) 

dz  X  - - - 

c  A 

U  =  cotj 


cosh  ( — ) 


1  ^max 

1  r-i  * 


Zsign(z') 


^  —  7 

^  “  ^min 


c(t-A)  _  -I  ,  U  +  cA.  ,  ,  u  x 

')  Y  [E(z',  t  -  u  / c)  -  E(z',  t  -  u  / c  -  A)  cosh"  ( - )-cosh  ( — ) 

]_  a,  a, 


u  =  ca 


It  is  important  to  note  that  the  summations  in  Ji(z,t)  and  J2(z,t)  involve  only  J(z',t’))  terms 
occurring  at  t'  <  t,  because  the  observation  point  is  separated  from  the  source  points  by  at  least 
the  gap  spacing  a. 

Evaluation  of  J3(z,t)  proceeds  in  a  manner  similar  to  J2(z,t)  but  with  a,  replaced  by  a0. 
Special  care  is  needed,  however,  to  treat  the  integrable  singularity  on  the  “self-patch”  where  z'  • 
z  and  f  « t.  Thus  we  write 

J3(z,t)  =  Js3(z,t)  +  J3ns(z,t) 
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where  the  "non-self'  term  follow  from  J2(z,t) 


1  Zmax  r  ,  .  .  .  .  ..1 

J“(z,t)  =  -  Esign(z')  E  [E(z  ,t-u/c)-E(z  ,t-u/c- A)Jx 


c(t-A) 


^  z'  =  z 


u  =  ca0  >  0 


cosh-1  (U  +  cA-)  -  cosh-1  (— ) 


a. 


a 


i  J 


a0  s  z-z 


and  the  self-patch  term  is  given  by 


d2.  rd,  ' 

7i  c<9t  Jz-cA/22  «p-z|  ^U2  -  (z'  - 

[E(z,l)-E(z,t-A)]2  f*'2  ft  1 
CA  7t  1  1  -y/u^ 


z)2 


=  [E(z,.)-E(z,1-A)]  2  ^ h_, (JL) 

cA  7i  J)  cA 

=  [E(z,t)-E(z,t-A)j  2  r  sech., (i)  +  ci  sjn-, (JL)' “ 
cA  re  L  cA  c A  Jo 

=  [E(z,  t)  -  E(z,  t  -  A)] 


“icA/2 


ln(2  +  V3)  |  1 
7t  3 


J(z,t)  J(M  A)  _  J[Em  (z  t)  _  Em  (Zjt  _  A)] 
V2t^3A 


m=l 


K, 


K  = 


ln(2  +  V3)  1 

+  3 


71 


Notice  that  the  self-patch  term  now  contains  a  term  in  J(z,t)  that  must  be  moved  to  the  left 
hand  side  of  the  MFIE.  Thus,  the  time-stepping  solution  can  be  written 


J inc  (z,  t)  +  J ,  (z,  t)  -  J2  (z,  t)-Jf  (z,  t)  +  K 


J(z,t)  = 


i^LA)+2Xtet)-E.(z,t-A)] 

V27tpA  ^ 


1  + 


K 


where  all  the  terms  on  the  right  involve  J(z't')  for  earlier  times  t'  <  t. 
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SECTION  6 


COMPARISON  WITH  ANALYTIC  SOLUTION 


In  this  section  we  will  validate  the  efficient  MFIE  solution  of  the  previous  section  against 
an  analytical  solution  of  the  same  problem.  We  find  good  agreement  bewteen  the  solutions  in 
both  time  and  frequency  domains. 


Ramo  Whinnery  and  van  Duzer  [20]  present  an  accessible  treatment  of  the  lossy  2D 
planar  transmission  line  in  the  frequency  domain.  For  a  wave  propagating  in  the  +z  direction  as 
exp(jcot-yz)  the  in  a  line  with  pc  =  p  and  sc  =  s,  the  propagation  constant  is  given  by 


i-j(j  +  D- 


where  8  is  the  skin  depth  at  frequency  co.  The  surface  impedance  boundary  condition  (SIBC) 
already  requires  5/a  «1,  so  the  propagation  constant  can  be  simplified  to 


Y  *  Jk 


i-Xj  +  i) 


8_ 

2a 


Thus,  an  initial  frequency  component  J0(oo)  e^01  becomes  J0(co)ej“‘  n  after  propagating  a  distance 

z  in  the  lossy  line.  The  real  part  of  y  causes  attenuation,  while  the  imaginary  part  gives  an 
increase  in  the  wave  number  that  reduces  the  propagation  speed.  We  then  use  the  inverse  fast 
Fourier  transform  to  obtain  the  analytic  solution  in  the  time  domain  for  comparison  to  the 
efficient  MFIE  solution.  Conversely,  we  also  performed  a  fast  Fourier  transform  of  the  efficient 
MFIE  solution  for  comparison  with  the  analytic  solution  in  the  frequency-domain. 


Before  proceeding  with  the  comparison  of  the  analytic  and  efficient  MFIE  solutions,  it  is 
important  to  relate  the  time-domain  good  conductor  requirement 

PA  » 1 

to  the  SIBC  (small  skin  depth  requirement)  and  another  requirement  used  in  the  derivation  of  the 
analytic  result  for  the  propagation  constant  y,  i.e., 


-  NN  1  CU1VJ  , - 

a  6  aV2 


respectively.  For  a  time-domain  computation  running  to  NtA,  the  relevant  frequency  range  is 


1 

NtA 


<f  <  — 
4A 


or 


2tc 

N^A 


<<£>< - 

tiA 


so  the  SIBC  and  analytic  requirements  respectively  become 


and 


VpA  » 


_a _ 1_ 

cA  n312 
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For  Nt  =  160  and  a  =  3cA  we  have  chosen  pA  =  1000  so  that  both  inequalities  are  satisfied  by 
least  a  factor  of  20.  It  is  important  to  note  that  the  efficient  integral  equation  approach  in  general 
requires  only  the  easily-satisfied  requirement 


VPA »-  1 


where  a  is  the  smallest  geometric  scale  length  of  the  problem. 


For  the  comparisons  we  have  redefined  the  physical  units  so  that  all  computational 
variables  are  dimensionless,  i.e.,  c=l,s=l,p=l,A=l.  The  driving  current  was  taken  to  be  a 
gaussian  modified  by  a  (1  +  cos)  factor  to  give  a  cutoff  at  t  =  40A, 


0.5  exp 


f  ( t  —  20A>\ 


2Y 


( 


1  +  COS 


j0(t)  =  j  t  8A 

[0  ,  otherwise. 

We  also  used  a  bipolar  shape  for  the  driving  current 


Vi 


71 


V 


20A 


—  1 


J) 


0  <  t  <  40A 


exp 


f  /t-20AV^' 


8A  ) 

\ 

0  ,  otherwise, 


f 


1  +  COS 


71 


V 


20A 


sm 


Y\ 


71 


20A 


-1 


,  0  <  t  <  40A 


J) 


which  was  then  normalized  to  give  a  peak  value  of  1 .  Although  both  the  monopolar  and  bipolar 
pulse  shapes  are  continuous  and  differentiable  at  t  =  0  and  t  =  40A,  the  higher  order  derivatives 
are  discontinuous  at  those  times. 


We  first  calculate  the  evolution  of  the  monopolar  pulse  in  an  ideal  (non-lossy)  line. 

Figure  6-1  shows  the  evolution  of  the  pulse  in  the  time  domain  as  it  passes  various  locations  in 
the  line.  As  expected  for  an  ideal  line,  the  pulse  propagates  with  no  change  in  shape.  Figure  6-2 
compares  the  magnitude  of  the  initial  frequency  spectrum  at  z  =  0  with  that  off  the  pulse  shape  at 
z  =  120cA.  We  find  excellent  agreement  out  to  a  frequency  of  about  0.16/A,  where  both  the 
initial  and  propagated  spectra  depart  from  the  smooth  decay  of  the  lower- frequency  spectrum. 
This  feature  is  due  to  the  temporal  cutoff,  which  as  we  noted  above  is  not  perfectly  smooth.  The 
transition  occurs  not  at  a  given  frequency,  but  rather  when  the  spectral  amplitude  reaches  a  low 
level  comparable  to  the  magnitude  of  the  time-domain  points  adjacent  to  the  cutoff. 
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Figure  6-1 .  Time-domain  propagation  of  a  monopolar  pulse  in  an  ideal  transmission  line  at 
locations  z  =  0,  40cA,  80cA,  &  120cA  using  the  MFIE  solution. 


Normalized  frequency,  fA 


Figure  6-2.  Frequency  spectra  of  a  monopolar  pulse  propagated  in  an  ideal  transmission  line  at 
z  =  0  (solid  line)  and  z  =  120cA  (dashed  line). 
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Next  we  treat  the  lossy  transmission  line,  for  which  we  use  the  bipolar  pulse.  The 
monopolar  pulse  has  a  zero-frequency  component,  which  has  an  infinite  skin  depth  that  violates 
the  SIBC  requirement.  We  have  chosen  pA  =  1000  and  10  terms  for  the  recursive  convolution. 
Figure  6-3  shows  good  agreement  between  the  analytic  and  efficient  MFIE  pulse  shapes  at 
various  positions  along  the  line.  As  expected  the  pulse  width  increases,  the  pulse  height 
decreases,  and  the  propagation  speed  (see  the  zero-crossing)  is  slightly  less  than  c  =  1 .  It  is 
important  to  note  that  both  analytic  and  efficient  MFIE  solutions  have  a  slowly  decaying  wake, 
so  that  any  finite-time  calculation  will  be  non-zero  at  t  =  NtA.  The  cutoff  of  this  wake  by  the 
finite  fast  Fourier  transform  will  always  introduce  artifacts  into  the  frequency  spectrum 

For  the  frequency  domain,  Figure  6-4  compares  the  analytic  and  efficient  MFIE  spectral 
amplitudes  at  z  =  0,  40A,  &  80A.  Again  the  agreement  is  good  at  lower  frequencies,  but  the 
temporal  cutoffs  cause  a  raised  spectral  plateau  at  high  frequency.  The  plateau  is  higher  for  z  = 
80cA  than  for  40cA,  because  the  former  waveform  has  less  time  than  the  latter  for  the  wake  to 
decay.  The  frequency-domain  attenuation  curves  provide  a  more  sensitive  comparison  of  the 
analytic  and  efficient  MFIE  solutions.  Figure  6-5  shows  good  agreement  at  frequencies  below 
the  cutoff,  which  is  really  an  artifact  of  the  finite  time  window  for  the  MFIE  computation. 

Thus  we  conclude  that  the  surface  impedance  boundary  condition  with  a  recursive 
convolution  is  a  valid  approach  for  the  solution  of  dispersive  time-domain  integral  equations. 


Figure  6-3.  Time-domain  propagation  of  a  bipolar  pulse  in  a  lossy  transmission  line  at  locations 
z  =  0, 40  cA,  80  cA,  &  120  cA;  solid  line  is  analytic  and  dashed  line  is  efficient  MFIE. 
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Figure  6-4.  Frequency  spectra  of  a  bipolar  pulse  propagated  in  a  lossy  transmission  line  at 
locations  z  =  0, 40  cA,  &  80  cA;  solid  line  is  analytic  and  dashed  line  is  efficient  MFIE. 


Normalized  frequency,  fA 


Figure  6-5.  Attenuation  of  a  bipolar  pulse  propagated  in  a  lossy  transmission  line  at  locations  z  — 
0,  40  cA,  &  80  cA;  solid  line  is  analytic  and  dashed  line  is  efficient  MFIE. 
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SECTION  7 


COMPUTATIONAL  CODES 


Several  computational  codes  have  been  developed  during  the  present  effort.  They  are  all  written 
in  Interactive  Data  Language  (IDL),  a  product  of  Research  Systems  Inc.  in  Boulder,  CO.  These 
codes  were  written  for  research  purposes  only,  and  are  not  intended  to  be  user-friendly  tools. 

LOSSYLINE  is  a  main  program  that  calculates  the  efficient  solution  of  the  dispersive  time- 
domain  integral  equation  for  a  lossy  2D  planar  transmission  line. 

;  main  program  lossyline 

n=160 

a=5. 

gw=8. 

pi=double(4.*atan(l .)) 
self=double(alog(2  .+sqrt(3  ,))/pi+ 1  ,/3 .) 
beta=double(1000.);  =timestep*sigma/epsilon/2 
conO=double(l  ,/sqrt(2.*pi*beta)) 

c=[1.2024728d0,2.5876867d-l,5.6953491d-2,1.2455136d-2,2.7406290d-3,  $ 
5.9659851d-4,1.3367602d-4,2.6443598d-5,8.4750592d-6,9.2249234d-7] 
w=[3  .dO,  1 .0908399d0,3 .96643  88d- 1 , 1 .4422497d- 1 ,5 .2442 1 1 5d-2,  $ 

1.9068649d-2, 6.9336146d-3, 2.521 1542d-3,9.1672520d-4,3.3333345d-4] 

mm=n_elements(w) 

cow=c/w/sqrt(beta) 

expw=exp(-w) 

jj=dblarr(n,n) 

con=dblarr(n,n) 

conm=dblarr(n,n,mm) 

incident, n, a, jj, gw 

for  t=l,n-l  do  begin 

print, t 

forz^Rt-l  do  begin 

if  (t  gt  1)  then  for  m=0,mm-l  do  conm(z,t,m)=conm(z,t-l,m)*expw(m)  $ 

+co  w(m)  *  exp  w(m)*  ( 1 .  -expw(m))  *jj  (z,t- 1 ) 
con(z,t)=total(float(conm(z,t,*))) 

jj(z,t)=(jj(z,t)+uzint(n,jj,a,z,t,con,conO)+uzint(n,jj,0.,z,t,con,conO)  $ 

+(j  j  (z,t- 1  )*  conO-(con(z,t)-con(z,t- 1 )))  *  self)/(  1  .+conO  *  self) 
endfor  ;  z 
endfor ; t 
ymin=-l 

plot ,j j (0, * ) ,yrange= [ymin,  1  ] , title-propagated  wave  vs  time' 

for  z=10,n-l,10  do  oplot,  jj(z,*) 

oplot,[0,n-l],[l,l] 
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end 


INCIDENT  is  a  subroutine  that  calculates  an  incident  waveform  for  either  the  LOSSYLINE  or 
IDEALLINE. 

pro  incident, n, a,jj,gw 
t=doubl  e(findgen(n)) 
pi=double(4.*atan(l .)) 

j  0=double(exp(-(t/gw-2 . 5)A2)*  ( 1  .+cos(pi  *(t/2. 5/  gw- 1 .)))) 

j  0=j  0  *  double(sin(-pi/2 .  *  (t/2 . 5  /  gw- 1 . ))) 

j0(0)=0. 

j0(5*fix(gw):*)=0. 

j0=2.*j0/maxG0) 

asq=a*a 

for  z=0,n-l  do  begin 
zsq=float(z)A2 
zal  =fix(sqrt(zsq+asq))+ 1 
fort=z,n-l  do  jj(z,t)=.5*h(t-z) 
for  t=zal  ,n-l  do  begin 
dum=double(0.) 
for  ul=zal  ,t  do  begin 
usq=float(ul)A2 
yl=usq-zsq 

if  (ul  eq  zal)  then  yO=asq  else  yO=(ul-l .)A2-zsq 
arg  1  =double((zsq-asq-2.  *  zsq*  asq/y  1 )/ (zsq+asq)) 
arg0=double((zsq-asq-2.*zsq*asq/y0)/(zsq+asq)) 

dum=dum+double(.  5  *  (( 1 .  -u  1 )  *j  0(t-u  1  )+u  1  *  j  0(t-u  1 + 1 ))  *  (asin(arg  1  )-asin(argO))  $ 

+z*  (j  0(t-u  1  )-j  0(t-u  1 + 1  ))*  (atan(sqrt(y  1  /asq- 1  .))-atan(sqrt(yO/asq- 1 .)))) 
endfor ;  ul 
jj(z,t)=jj(z,t)-dum/pi 
endfor ;  t 
endfor ;  z 
return 
end 

UZINT  is  a  function  that  calculates  the  dz’  du’  integral  for  LOSSYLINE. 

function  uzint,n,jj, a, z,t, con, conO 
zzmax=fix((float(t)A2-float(z)A2-aA2)/2./(t-z)) 
zzmin=(-fix((float(t)A2-float(z)A2-aA2)/2./(t+z))) 
duml=0. 

for  zz=zzmin,zzmax  do  begin 
zzz=abs(zz) 

alpha=double(sqrt(aA2+float(zz-z)A2)) 

dum=0. 

cl=0. 
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dl=0. 

for  uO=fix(alpha),t-l  do  begin 

ul=float(uO)+l. 

cO=cl 

dO=dl 

c  1  =double(sqrt(u  1  A2-alphaA2)) 
if  (alpha  le  0.)  then  dl=0.  else  $ 

dl=double(alog(ul/alpha+sqrt((ul/alpha)A2-l.))) ;  asech(alpha/ul) 

tau=t-u0 

jjj=jj(zzz,*) 

;  lossy  convolution  contribution 
if  (uO  ge  1)  then  dum-dum-  $ 

(((jjj  (tau)-jjj  (tau- 1  ))*conO-(con(zzz.tau)-con(zzz,tau- 1  )))*  (d  1  -d0)) 

;  non-lossy  contribution 

if  (a  ge  1.)  then  dum=dum+double(a*(jjj(tau)-jjj(tau-l))*(cl-cO)/alphaA2) 
endfor ;  u 

duml:=duml+dum*(fix(zz  gt  0)-fix(zz  It  0)) 
endfor ;  zz 

return,  duml/double(4.*atan(l .)) 
end 

IDEALLINE  is  a  main  program  that  calculates  a  a  time-domain  integral  equation  solution  for  the 
ideal  2D  planar  transmission  line. 

;  main  program  idealline 

n=160 

gw=8. 

a=2. 

jj=fltarr(n,n) 
incidents, a, jj, gw 
for  t=l  ,n-l  do  begin 
print, t 

for  z=  1  ,t- 1  do  j  j  (z,t)=jj  (z,t)+uzinti(n,jj  ,a,z,t) 

endfor ; t 

ymin=0 

plot, jj(*,0),yrange=[ymin,l], title-propagated  wave  vs  distance' 
for  t=10,n-l,10  do  oplot, jj(*,t) 
oplot,  [0,n- 1],[1,1] 

plot, jj(0,*),yrange=[ymin,l], title-propagated  wave  vs  time' 
for  z=10,n-l,10  do  oplot,  jj(z,*) 
oplot,  [0,n-l],  [1,1] 
end 


UZINTI  is  a  function  that  returns  the  dz’  du’  integral  for  LOSSYLINE. 
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function  uzinti,njj,a,z,t 

zzmax=fix((float(t)A2-float(z)A2-aA2)/2 ./  (t-z)) 

zzmin=(-fix((float(t)A2-float(z)A2-aA2)/2./(t+z))) 

duml=0. 

for  zz=zzmin,zzmax  do  begin 
zzz=abs(zz) 

alpha=sqrt(aA2+float(zz-z)A2) 

dum=0. 

cl=0. 

for  uO=fix(alpha),t-l  do  begin 

ul=float(uO)+l. 

cO=cl 

c  1  =sqrt(u  1  A2-alphaA2) 
tau=t-uO 

jjj=jj(zzz»*) 

dum=dum+a*  (jjj  (tau)-jjj  (tau- 1  ))*  (c  1  -cO)/alphaA2 
endfor ;  uO 

duml=duml+dum*(fix(zz  gt  0)-fix(zz  It  0)) 
endfor ;  zz 

return,  dum  1/(4.* atan(  1 .)) 
end 

ANALYTIC  is  a  main  program  that  compares  the  analytic  and  computational  solutions.  It  runs 
immediately  after  LOSSYLINE  and  assumes  that  the  variables  jj  and  beta  area  available. 

;  main  program  analytic 

pi=4.*atan(l.) 

dt=l. 

t=dt*findgen(n) 

freq=indgen(n) 

n21=n/2+l 

freq(n2 1  )=n2 1  -n+findgen(n2 1  -2) 

freq=freq/float(n)/dt 

w=2.*pi*freq 

f=freq(0:n/2) 

i=complex(0.,l.) 

k=w  ;  wavenumber  for  speed  of  light  =  1 . 
u=[0.,sqrt(2.*abs(w(l  :*))/(2.*beta))/a/abs(k(l  :*))] 
aa=sqrt(sqrt(  1  .+2.*u+2.*uA2)+l  +u) 
gamma=sqrt(abs(w)/(2 .  *beta))/ a/aa+i  *  k*  aa/sqrt(2 .) 
j0=jj(0,*) 
specO=fft(jO,-l) 

spec40a=spec0  *  exp(-4  0  .*  gamma) 
spec  80a=spec0  *  exp(-  8  0 .  *  gamma) 
sped  20a=spec0*exp(-l  20.*gamma) 
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j40a=float(fft(spec40a,  1 )) 
j  80a=float(fft(spec80a,  1 )) 
j  1 20a=float(fft(spec  1 20a,  1 )) 
spec40c=fft(jj  (40,*),- 1 ) 
spec80c=fft(ij(80,*),-l) 

plot_io,f,abs(spec0),xrange=[0,.20],yrange=[le-7,.l],yticks=6,  $ 

title='abs(spectrum)  vs  frequency' 

oplot,f,abs(spec40a) 

oplot,f,abs(spec80a) 

oplot,f,abs(spec40c),linestyle=9 

oplot,f,abs(spec80c),linestyle=9 

plot, t,j0,xrange=[0, 160.], xticks=8,yrange=[-l,l],$ 

title='Propagated  wave  vs  time  at  z=0,40,80,120' 

oplot,t,j40a 

oplot,t,j80a 

oplot,t,jl20a 

oplot,t,jj(40,*),linestyle=9 
oplot,t  jj  (80,  *  ),linestyle=9 
oplot,t,jj(l  20,*),linestyle=9 
db=20./alog(10.) 

plot,f,db*80.*float(gamma),xrange=[0,.2],  $ 

title='Attenuation  (dB)  vs  frequency  at  z=40,80' 

oplot,f,db*40.*float(gamma) 

oplot,f,db*alog(abs(spec0/spec40c)),linestyle=9 

oplot,f,db*alog(abs(spec0/spec80c)),linestyle:=9 

end 

EXPFIT  is  a  main  program  that  fits  Zi  with  a  sum  of  exponentials. 

;  main  program  expfit 

pi=4.*atan(l.) 

n=100 

m=10 

q=3. 

tmax=l.e3 

tmin=l. 

time=tmin*exp(alog(tmax/tmin)*fmdgen(n)/float(n-l)) 

zz=l  ,/sqrt(pi*(2.*time)A3)  ;*(1  ,+3./8/time+45./128./timeA2+525./1024./timeA3) 

w=double(q/tmin*exp(-findgen(m)*alog(tmax/tmin*qA2)/float(m-l))) 

expw=exp(-w) 

a=dblarr(m,m) 

b=dblarr(m) 

f=dblarr(m,n) 

for  i=0,m-l  do  begin 

for  k=0,n-l  do  f(i,k)=exp(-w(i)*time(k))/zz(k) 
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b(i)=total(f(i,*)) 

endfor 

for  i=0,m-l  do  begin 

for  j=0,m-l  do  a(i,j)=total(f(i,*)*f(j,*)) 

endfor 

bb=max(b) 

c=cramer(a/bb,b/bb, /double) ;  an  IDL  library  routine  to  solve  linear  equations 
for  i=0,m-l  do  print, i,c(i),w(i),l./w(i) 
print, total(c),total(c*w),total(c/w) 
z=dblarr(n) 

for  k=0,n-l  do  z(k)=total(c(*)*exp(-w(*)*time(k))) 
plot_oo,time,zz,title='Impedance  Z1  vs  time:  exact  and  fit' 
oplot,  time,  z,  linestyle=9 

plot_oi, time, z/zz,yrange=[. 99,1 .01], title='(exP  fit)/ exact  relative  error  vs  time' 
oplot,  [tmin,tmax],[  1 , 1  ] 
print, max(zz/z),min(zz/z) 
end 
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